clear all
t_oborota = 12*3600;%vremya oborota v sekyndah
c = 3e5;%skorost sveta
f_0 = 1.602e9;%chastota glonass
T = 60; %shag modeli sec
R_zem = 6.4e3;%radius zemli
R_orb = 19.1e3;%visota zemli
R_sv = R_orb + R_zem;%sputnik
t = 0:T:360*T;
alfa = (2*pi/t_oborota)*t;
Rsv_x = R_sv * cos(alfa); 
Rsv_y = R_sv * sin(alfa);
subplot(2,1,1);
plot(t,Rsv_x,t,Rsv_y);
subplot(2,1,2);
plot(Rsv_x,Rsv_y);
Rr = [0,R_zem];
%vizirovanie
Rvz_x = Rsv_x - Rr(1);
Rvz_y = Rsv_y - Rr(2);
figure(2);
plot(t,Rvz_x,t,Rvz_y);
%skorosti sputnika
V_x = diff(Rsv_x)/T;
V_y = diff(Rsv_y)/T;
V = sqrt(V_x.^2 + V_y.^2);
R_abs = sqrt(Rvz_x.^2 + Rvz_y.^2);
figure(3);
plot(V_x,V_y);
%proekciya skorosti
V_abs_vz = -(V_x.*Rvz_x(2:end)+V_y.*Rvz_y(2:end))./(R_abs(2:end));
%vi4islenie doplera
f_d = V_abs_vz.*f_0/c;
figure(4);
plot(f_d);
%massiv vidimosti Rvz_y>0
%massiv nevidimosti Rvz_y<0
k_1 = 0; k_2 = 0;
for k= 1:(length(t)-1)
    if Rvz_y(k)>0
        k_1 = k_1+1;
        f_d_v(k_1) = f_d(k);
        x_v(k_1) = Rsv_x(k);
        y_v(k_1) = Rsv_y(k);
    else 
        k_2 = k_2+1;
        f_d_nv(k_2) = f_d(k);
        x_nv(k_2) = Rsv_x(k);
        y_nv(k_2) = Rsv_y(k);
    end;
end;
figure(5);
plot(x_v,y_v,x_nv,y_nv);
figure(6);
m=-5e3:100:5e3;
hist(f_d_v,m);


